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ABSTRACT 

Ill  this  paper  we  apply  a  sensitivity  equation  method  to  shape  optimization  problems. 
An  algorithm  is  developed  and  tested  on  a  problem  of  designing  optimal  forebody  simulators 
for  a  2D,  inviscid  supersonic  flow.  The  algorithm  uses  a  BFGS/Trust  Region  optimization 
scheme  with  sensitivities  computed  by  numerically  approximating  the  linear  partial  differ¬ 
ential  equations  that  determine  the  flow  sensitivities.  Numerical  examples  are  presented  to 
illustrate  the  method. 
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1  Introduction 


The  development  of  practical  computational  metho<ls  for  optimization  based  design  and 
control  often  relies  on  cascading  simulation  software  into  optimization  algorithms,  [^lack-box 
methods  are  examples  of  this  approach.  Although  the  precise  form  of  the  overall  “optimal 
design”  (OD)  algorithm  may  change,  there  is  an  often  unstated  assumption  that  properly 
combining  the  “best”  simulation  algorithm  with  the  “best”  optimization  scheme  will  produce 
a  good  OD  algorithm.  There  are  many  examples  to  show  that  in  general  this  assumption  is 
not  valid.  However,  in  many  cases  it  is  a  valid  assumption  and  often  this  approach  is  the 
only  practical  way  of  attacking  complex  optimal  design  problems.  If  one  uses  this  cascading 
approach,  then  it  is  still  important  to  carefully  pass  information  between  the  simulation 
and  the  optimizer.  Typically,  one  uses  a  simulation  code  to  produce  a  finite  dimensional 
model  and  this  discrete  model  is  then  used  to  supply  approximate  function  evaluations  to 
the  optimization  algorithm.  Moreover,  the  approximate  function.'^  are  then  differentiated  to 
supply  gradients  needed  by  the  optimizer.  Although  there  are  numerotis  variations  on  this 
theme,  they  all  may  be  formulated  as  “approximate-then-optimize”  approaches.  There  are 
other  approaches  that  first  formulate  the  problem  as  an  infinite  dimensional  optimization 
problem  and  then  use  numerical  schemes  to  approximate  the  optimal  design.  All-at-once, 
one-shot  and  adjoint  methods  are  examples  of  this  “optimize-then-approximate”  approach. 
Regardless  of  which  approach  one  chooses,  some  type  of  approximation  must  be  introduced 
at  some  point  in  the  design  process. 

The  sensitivity  equation  (SE)  method  is  an  approach  that  views  the  simulation  scheme 
as  a  device  to  produce  approximations  of  both  the  function  and  the  sensitivities.  The  ba¬ 
sic  idea  is  to  produce  approximations  of  the  infinite  dimensional  sensitivities  and  to  pass 
these  “approximate  derivatives”  to  the  optimizer  along  with  the  approximate  function  eval¬ 
uations.  There  are  several  theoretical  and  practical  issues  that  need  to  be  considered  when 
this  approach  is  u.sed.  For  example,  there  is  no  assurance  that  the  SE  method  produces  “con¬ 
sistent  derivatives.”  This  will  depend  on  the  particular  numerical  scheme  used  to  discretize 
the  problem.  However,  the  SE  method  allows  one  the  0])tion  of  using  separate  numerical 
schemes  for  flow  solves  and  sensitivities,  so  that  consistent  derivatives  can  be  forced.  VVe 
shall  not  address  these  issues  in  this  short  paper.  The  goal  here  is  to  illustrate  that  a 
SE  based  method  can  be  used  with  stamlard  optimization  schemes  to  produce  a  practical 
fast  algorithm  for  optimal  design.  We  concentrate  on  a  particular  application  (the  optimal 
forebody  design  problem)  and  u.se  a  specific  iterative  solver  for  the  flow  equations  (F'AR('). 
Many  flow  solvers  are  itc’iative  and  for  the.se  types  of  codes,  the  SE  method  has  perhaps  the 
maximum  potential  for  improving  speed  and  accuracy. 


In  the  next  section  we  describe  tlie  forebody  design  problem  and  formulate  the  oiitimal 
design  problem.  In  Sections  3  and  4  we  review  the  derivation  of  the  sensitivity  equations 
and  in  Section  5  we  discuss  modifications  to  an  existing  simulation  code  that  are  needed  in 
order  to  use  that  code  for  computing  sensitivities.  In  Section  6,  we  present  numerical  results 
for  the  optimal  design  problem  and  Section  7  contains  conclusions  and  suggtvstions  for  intuit' 
work. 


2  Optimal  Design  of  a  Forebody  Simulator 

This  problem  is  a  2D  version  of  the  problem  described  in  [1,4,8]  .  The  Arnold  Engineering 
Development  fVuter  (AED(l)  is  developing  a  free-jet  test  tacility  for  full-scale  testing  of 
engines  in  various  free  flight  conditions.  Although  the  test  cells  are  large  enough  to  house 
the  jet  engines,  they  are  too  small  to  contain  the  full  airplane  forebody  and  engine.  Thus, 
the  effect  of  the  forward  fuselage  on  the  engine  inlet  flow  conditions  must  be  “simulated.” 
One  approach  to  solving  this  problem  is  to  replace  the  actual  forebody  by  a  smaller  object, 
called  a  “forebody  simulator”  (FBS),  and  determine  the  shape  of  the  FB.S  that  produces  the 
best  flow  match  at  the  engine  inlet.  The  2D  version  of  this  problem  is  illustrated  in  Figure 
2.1  (see  [l],[4],  [8]  and  [9]). 

The  underlying  mathematical  model  is  based  on  conservation  laws  for  mass,  momentum 
and  energy.  For  inviscid  flow,  we  have  that 
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The  velocity  components  u  and  n,  the  pressure  F,  the  temperature  T,  and  the  Mach  number 
M  are  related  to  the  conservation  variables,  i.e.,  the  components  of  the  vector  Q,  by 
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At  the  inflow  boundary,  we  want  to  simulate  a  free-jet,  so  that  we  specify  the  total 
pressure  Fq,  the  total  temperature  Tq  and  the  Mach  number  Mq.  We  also  set  7>  =  0  at  the 


2 


inflow  bouiulary.  If  tt/,  Pj  and  Ty  denote  the  inflow  values  of  the  x-coinponent  of  the  velocity, 
the  pressure  and  the  temperature,  these  may  be  recovered  from  Po,  To  and  Mo  by 


r,  = 


To 


Pi  = 


Po 


and  uj  =  MqTj  = 


M^To 


(4) 

The  components  of  Q  at  the  inflow  may  then  be  determined  from  (4)  through  the  relations 


'1'^/  n  J  r  Pi  ,  trs 

Pi  =  -TjfT-,  fill  =  PiUh  ni  =  0  and  Ej  = - -  +  .  (5) 
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The  forebody  is  a  solid  surface,  so  that  the  normal  component  of  the  velocity  vanishes, 

i.e., 

uj/i  +  vt)2  =  0  on  the  forebody,  (6) 

where  7;,  and  7/2  are  the  components  of  the  unit  normal  vector  to  the  boundary.  Note  that 
we  impose  (6)  on  the  velocity  components  u  and  v,  and  not  on  the  momentum  components 
777  and  77.  Insofar  as  the  state  is  concerned,  it  is  clear  that  it  does  not  make  any  difference 
whether  (6)  is  imposed  on  m  and  n  or  on  u  and  v,  since  m  =  pu  and  n  =  pv  and  p  ^  0.  It 
can  be  shown  that  it  does  not  make  any  difference  to  the  sensitivities  as  well. 

Assume  that  at  x  =  the  desired  steady  state  flow  Q  =  Q{y)  is  given  as  data  on  the  line 
(called  the  Inlet  Reference  Plane) 


IRP  =  {(x,y)|x  =  0,(T<y  <  6}. 


Also,  we  assume  here  that  the  inflow  (total)  Mach  number  Mo  can  be  used  a.s  a  design 
(control)  variable  along  with  the  shape  of  the  forebody.  Let  the  forebody  be  determined  by 
the  curve  P  =  r(x),  a  <  x  <  ft  and  let  p  =  (Mo,r(-)).  The  problem  can  be  stated  as  the 
following  optimization  problem: 

Problem  FBS.  Cl  iven  data  Q  =  Q{y)  on  the  I RP,  find  the  parameters  p*  =  (A/q,  r*(-)) 
such  that  the  functional 

=  \\Qoo{ft,y)-Q{y)fdy  (7) 

is  minimized,  where  Qoo{Xi  y)  =  Qoo{x,  y,  p)  is  the  solution  to  the  steady  state  Euler  equations 

In  the  FB.S  design  problem,  the  data  Q  is  generated  both  experimentally  and  numerically. 
In  particular,  the  full  airplane  forebody  (which  is  longer  and  larger  than  the  desired  FBS)  is 
used  to  generate  the  data.  Since  the  FBS  is  “constrained”  to  be  shorter  and  smaller,  we  shall 


consider  the  optimization  problem  illustrated  in  Figure  2.2  below.  The  data  Q  is  generated 
by  solving  (l)-(6)  for  the  long  forebody  in  Figure  2.2-(a)  and  the  problem  is  to  find  p’  to 
minimize  where  the  shortened  FBS  is  constrained  to  be  one  half  the  length  of  the  “real 
forebody.”  This  problem  provides  a  realistic  test  of  the  optimal  design  algorithm  in  that  the 
data  can  not  be  fitted  exactly.  Also,  we  note  that  we  have  a  problem  with  shocks  in  the  flow 
field.  As  shown  in  [2],  optimization  of  flows  with  shocks  can  be  difficult  and  requires  some 
understanding  of  the  impact  that  shocks  have  on  the  smoothness  of  the  cost  functional. 

(dearly  the  statement  of  the  problem  is  not  complete.  For  example,  one  should  carefully 
specify  the  set  of  admissible  curves  r(  )  and  questions  remain  about  existence,  uniqueness 
and  integrability  of  “the”  solution  Qoo-  We  will  not  address  these  issues  in  this  short  paper. 

Most  optimization  based  design  methods  require  the  computation  of  the  derivatives 
-^Q^{x,y,p).  These  derivatives  are  called  sensitivities  and  various  schemes  have  been  de¬ 
veloped  to  approximate  the  sensitivities  numerically  (see  [7],  [8],  [10]  and  [11]).  A  common 
approach  is  to  use  finite  ditferences.  In  particular,  the  steady  state  equation  (8)  is  solved 
for  p  and  again  for  p  +  Ap  and  then  -§^Q ooi-i: ,  y ,  P)  is  approximated  by  . 

This  method  is  often  costly  and  can  introduce  large  errors.  Another  approach  is  to  first 
derive  an  equation  (the  sensitivity  equation)  for  ooi^ ,  V ,  p)  and  then  numerically  solve 
this  equation.  We  shall  illustrate  this  approach  for  the  forebody  design  problem.  In  the  next 
two  sections  we  derive  the  sensitivity  equations.  Although  these  derivations  may  be  found 
in  [3]  we  repeat  them  here  for  completeness. 


Sensitivities  with  Respect  to  the  Inflow  Mach  Num¬ 
ber 


First,  we  consider  the  design  parameter  Mq.  We  will  derive  equations  for  the  sensitivity 
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(10) 


dM^'  ~  dMo^'  ~  dM^  ^  ~  dMl 

The  differential  equation  system  (1)  has  no  explicit  dependence  on  the  design  parameter 
Ml  so  that  equations  for  the  components  of  Q'  are  easily  determined  by  formally  differen¬ 
tiating  (1)  with  respect  to  Ml  The  result  is  the  system 


^  dF[ 

dt  dx  dy 
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(11) 
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where 


mu'  +  m'u  +  P' 
mv'  +  m'v 

(£  +  PK  +  (E'+P')u 


and  F2  = 


nu'  +  7i'u 
nv'  -i-  7i'v  +  P' 

{E  +  P)v'  +  {E'  +  P')v 


and  where, 


dMi' 


and  r  = 


dMi' 


and  where,  through  (3),  the  sensitivities  (10)  and  (13)  are  related  by 

u'  =  -771  -  —p,  P'  =  {'>(-  \)(e'  -  ^/)'(li’^  +  v^)  -  p{uu  +  vv'))  , 

p  p-‘  \  I  / 

v'  =  —71' - -p'  and  T'  =  7(7  —  1)  ( -E' - -p'  —  iuu'  +  vv')]  .  (14) 

p  p^  \p  p^  ) 

Note  that  (11)  is  of  the  same  form  a^i  (1),  with  a  different  flux  vector.  In  particular,  (11) 
is  in  conservation  form.  As  a  result  of  the  fact  that  (11)  is  Imear  in  the  primed  variables, 
and  that  by  (14)  n',  v'  and  P'  are  linear  in  the  components  of  Q' ,  (11)  is  a  linear  system  in 
the  sensitivity  (9),  i.e.,  in  the  components  of  Q' . 

Now,  we  need  to  discuss  the  boundary  conditions  for  Q'.  Except  for  the  inflow  conditions, 
all  boundary  conditions  are  independent  of  the  design  parameter  Mq.  Thus,  the  latter  may 
be  differentiated  with  respect  to  Mq  to  obtain  boundary  conditions  for  the  sensitivities.  For 
example,  at  the  forebody  where  (6)  holds,  we  simply  would  have  that 

u't]\  +  v'tj2  =  0  on  the  forebody.  (15) 

Similar  operations  yield  boundary  conditions  for  the  sensitivities  along  symmetry  lines,  other 
solid  surfaces  and  at  the  outflow  boundary.  Note  that  if  instead  of  (6),  one  interprets  the  no 
penetration  condition  as  one  on  the  momentum,  i.e.,  7717]^  +  7i7}2  =  0  on  the  forebody,  then 
instead  of  (15)  we  would  have  that 

m'r/i  +  717]2  =  0  on  the  forebody  (16) 


which  is  seemingly  different  from  (15).  However,  (6)  and  (14)  can  be  used  to  show  that 


m'r/i  +  7l'7f2  =  p{u'7]i  +  v'TJi)  +  p'{u7]^  +  V7]2)  =  p{u'7]i  +  (17) 

so  that,  since  p  ^  0,  (15)  and  (16)  are  identical. 

The  inflow  boundary  conditions  for  the  sensitivities  may  be  determined  by  differentiating 
(4)  and  (5)  with  respect  to  the  design  parameter  Mq.  Note  that  this  parameter  appears 


explicitly  in  the  right-hand-sides  of  the  eciuations  in  (4)  and  (5).  Without  difficulty,  one 
finds  from  (5)  that 
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where,  from  (4), 
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4  Sensitivities  with  Respect  to  the  Forebody  Design 
Parameters 


We  assume  that  the  forebody  is  described  in  terms  of  a  finite  number  of  design  parameters 
which  we  denote  by  Pk,  k  =  1, . . . ,  /v,  and  that  the  forebody  may  be  described  by  the  relation 


y  =  ^{x;PuP2,---,Pk),  a  <  x  < 


(20) 


We  express  the  dependence  of  the  state  variable  Q  on  the  coordinates  and  the  design 
parameters  by  Q  =  Q{t,  x,  y\  Mq,  P\,  Pi^. . .  Pk)-  We  have  already  seen  what  eciuations  can 
be  used  to  determine  the  sensitivity  of  the  state  with  respect  to  i.e.,  for  Q'.  We  now 
discuss  what  equations  can  be  used  to  determine  the  sensitivities  with  respect  to  the  forebody 
design  parameters  Pk,  k  =  ,  K,  i.e.,  for 
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System  (1)  has  no  explicit  dependence  on  the  design  parameters  Pk,  so  that  equations 
for  the  components  of  Qk  are  easily  determined  by  differentiating  (1)  with  respect  to  Pk, 
k  —  1, . . . ,  ii.  This  produces  the  systems,  k  =  1, . . . ,  A' ,  given  by 
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dt  dx  dy 
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where 
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and  where, 
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Moreover,  by  (.3),  the  sensitivities  (22)  and  (25)  are  related  by 

Uk  =  yUk  -  ypk.  Pk  =  (7  -  I  )  +  VVk)^  , 

1  7t  ( \  E  'N 

Vk  =  ~nk - and  T*,.  =  7(7  -  1)  ( -F/t - ::Pk  -  {uuk  +  Wk)  \  ,  (26) 

p  p^  \p  r  J 

for  k=l,. . .  ,K. 

All  boundary  conditions  except  the  one  on  the  forebody  also  do  not  depend  on  the 
forebody  design  parameters  Pk,  k  =  For  example,  consider  the  inflow  boundary 

conditions  (4)-(5).  Differentiating  these  with  respect  to  Pk,  k  =  I, . . . ,  K  yields  that 

pki  =  ”U-/  =  ^iki  =  Eki  =  Tkt  =  Pki  =  Uki  =  Vkt  =  0  (27) 


at  the  inflow  boundary.  Now  consider  the  boundary  condition  (6)  on  the  forebody.  We  have 
that  on  the  forebody 

^  (28) 

7/2  Ox 


(’ombining  (6)  and  (28)  we  have  that 


u— - u  =  0 

Ox 


along  the  forebody  or,  displaying  the  full  functional  dependence  on  the  coordinates  and 
design  parameters,  we  have  at  a  point  (x,y)  on  the  forebody,  and  at  any  time  t, 

U  {t,X,y  =  ^X-,P„P2,...,PH);MlPuP2,....PK)-^{x-,Px,P2,....Ps) 

-V  (t,  X,  y  =  <D(x;  P, ,  P2,  •  -  • ,  Pa  );  M'^,  P,  ,  P2, .  • . ,  Pa  )  =  0.  (30) 

We  can  proceed  to  differentiate  (30)  with  respect  to  any  of  the  forebody  design  parameters 
Pk,  k  =  1, . . . ,  A' .  The  result  is  that,  along  the  forebody  for  A;  =  1, . . . ,  A', 
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where  u,  e  aiul  their  derivatives  are  evaluated  at  the  forebody  (x,j/  = 

If  an  iterative  scheme  is  used  to  find  a  steady  state  solution  of  this  system  ((23),  (27), 
(31)),  then  we  assume  that  present  guesses  for  the  state  variables  u  and  v  and  their  deriva¬ 
tives  ()u/dy  and  dvjdy  and  for  the  design  parameters  Mq  and  k  —  1, .  . .  ,  A’,  are  known. 
It  follows  that  the  right-hand-side  of  (31)  is  known  as  well  and  equation  (31),  the  bound¬ 
ary  conditions  along  the  forebody  for  the  sensitivities  with  respect  to  the  forebody  design 
parameters,  is  merely  an  inhomogeneous  version  of  (29),  the  boundary  condition  along  the 
forebody  for  the  state. 

Let  us  now  specialize  to  the  type  of  forebodies  considered  by  Htiddleston,  [8,9],  i.e.. 


<D(x;P,,P„...,Pa)  =  ^P,<?,(x), 

*.=] 


(.32) 


wher  ,  .,(r),  k  =  1,...,  A',  are  prescribed  functions,  e.g.,  Bezier  curves  (see  [6]).  In  this 


case. 


and 


Jw, 


and 


(:{:?) 


(34) 


('ombining  (3I)-(34),  one  obtains  that,  at  any  point  (j:,4>(x))  on  the  forebody  and  for  eacli 
k  =  1,...,  A, 


n 

E  I  = 

u=l 


,j=i 


dx 


dx 


(35) 


For  forebodies  of  the  type  (32),  (35)  gives  the  boundary  conditions  along  the  forebody  for 
the  sensitivities  with  respect  to  the  forebody  design  parameters  P*,-,  A;  =  1, . . . ,  A'.  It  is  now 
clear  that,  given  guesses  for  the  state  v?*riables  u  and  v  and  their  derivatives  dti/dy  and 
dv/dy  and  for  the  design  parameters  Mq  and  P*..,  k  =  1,...,  A',  then  the  right-hand-side  of 
(35)  is  known. 

(’onsider  now  the  problem  of  minimizing  J{p)  as  defined  above.  Most  optimization 
algorithms  use  gradient  information.  In  particular,  if  P*.  denotes  one  of  the  shape  parameters, 
then  the  derivative 

=  l  <  [^Qoo(/i,J/,/5) 

may  be  required  in  the  optimization  loop.  The  sensitivity  g^Qoc>(^i  J/i  p)  sati.sfies  the  steady- 
state  version  of  the  sensitivity  equations  (23).  In  practice  one  must  construct  approximations 
p)  and  feed  this  information  into  the  optimizer. 


,Qoo{fi,y,p)  -  Q{y)  >  dy  (.36) 
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Assume  that  one  has  a  particular  simulation  srheine  (finite  differences,  finite  elements, 
etc.)  to  approximate  the  How  on  a  given  grid,  i.e. 


Qh{‘r^y,p)  Qoc.(j-,y,p). 


CM] 


as  the  “step  size”  h  0.  (Jiven  the  design  parameter  p,  one  constructs  a  grid  (depending 
on  p)  and  then  omputes  Qll{J‘^,y^p)  ~  This  process  may  recpiire  some  type  of 

iterative  sclu'ine.  We  will  address  this  issue  below.  In  theory,  one  could  use  the  same  grid  and 
computational  scheme  to  approximate  y,p)  so  that  one  generates  •‘a])proximate 

sensitivities” 


() 


■:^Qoc.{j:^y-P) 

at  k 


d 

i)Pk 


Q.^.(x,y,p) 


as  /j  — >  0.  It  is  important  to  note  that  in  general 


^Qoc,(-r,y,p)  ^  ^lQhi-r,y,p)], 


(;w) 


(:W) 


i.e.  this  approach  may  not  provide  “consistent  sensitivities”.  However,  some  schemes  do 
provide  consistent  derivatives  and  even  if  (39)  holds,  the  error 


EDn  = 


df\ 


Q^{x,y,p) 


d 


[Qu{j-.y,p)] 


(40) 


I,, 

may  be  sufficiently  small  so  that  the  optimization  algorithm  converges.  Trust  region  methods 
are  particularly  well  suited  for  problems  of  this  type,  where  derivative  information  may  con¬ 
tain  (small)  errors.  As  we  shall  see  below,  there  are  certain  cases  where  f'an 

be  computed  fast  and  accurately.  Hence,  the  SE  method  provides  estimates  for  sensitivitir's 
that  may  prove  “good  enough”  for  optimization  and  yet  relatively  cheap  to  compute.  A  com¬ 
parison  of  7j7^Qoc.(j'i  ,</,  p)]^  and  various  finite  difference  approximations  of  ^  [Q/,(x,  y,  p)) 
may  be  found  in  [3]. 

It  is  important  to  note  that  the  details  of  the  computations  needed  to  approximate  a 
sensitivity  are  not  the  central  issue  here.  For  example,  the  sensitivity  equations  (11)  and 
(23)  are  viewed  as  independent  partial  differential  equations  that  must  be  solved  by  “some” 
numerical  scheme.  This  scheme  does  not  necessarily  have  to  be  the  same  scheme  \ised  to 
solve  the  flow  equation  (1),  although  as  we  shall  see  below,  there  are  ca.ses  where  using  the 
same  scheme  is  a  useful  approach. 

Also,  note  that  the  .sensitivity  equations  are  derived  for  the  jiroblem  formulated  on  the 
“physical”  domain.  If  one  tises  a  computational  method  that  maps  the  problem  to  a  com¬ 
putational  domain  (as  does  PAR('),  then  the  SE  method  does  not  require  derivatives  of 
this  mapping.  One  simply  maps  the  sensitivity  equation  (including  the  necessary  bound¬ 
ary  conditions),  grids  the  computational  domain,  solves  the  resulting  transformed  equations 
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him!  tluMi  maps  l»ark  to  t  lif  physical  (iomaiii.  If.  ou  th**  other  liaii<l.  one  iiiappeil  I  lie  flov\ 
eipiatioii  ( 1 )  ainl  rleriverl  a  sensitivity  ecpiatioii  in  the  <'omi;ntational  domain,  then  to  obtain 
the  correct  sensitivities  one  would  have  tc»  compute  the  mapping  sensitivity.  1  her<  fore.  it  i> 
m<»re  efficient  tcj  derive  the  sensitivity  erpiatioiis  in  the  physii  al  domain. 

Finally,  we  not*’  that  th*’  SK  inetho*!  <l«’scrilie<|  here  has  on*’  a<iditi*)nal  hca  fit.  lo  com 
put*’  a  .s*’n.sitivity.  say  t/tQ,  {x.  ij.p).  tluni  on*’  first  s*’l«’cts  th*’  param«’t*’r  vain*-  /».  constnn  ts 
a  ( i>n;putati*jnal  grid  and  solv*‘s  f*jr  ,  (r.  »/.p)j^ .  l  ln-r*’  is  no  n*’«’d  to  (  (uiipute  gricl 

s*'nsitiviti«’s. 

5  Computing  Sensitivities  using  an  Existing  Code  for 
the  State 

Supp<j.s«’  on*’  ha,s  available  a  co*le  t(j  compute  th*’  stat*’  variable’s,  i.e..  to  fiinl  appnj.ximate 
s*>lutions  ( 1  )  ahuig  with  bouinlary  and  initial  < oinlitions.  In  principle,  it  is  an  *’asy  matter 
to  ann’iid  such  a  c(jde  so  that  it  can  al.s*j  comput**  s<’nsitiviti*‘s. 

First.  I*’t  us  <'ompare  (I )  with  (II).  If  <uie  wishes  t*>  atn<’n*l  th**  «’.xisting  xulr  that  r  an 
hainlle  (1)  so  that  it  can  treat  (11)  as  w<‘ll,  *jn*’  has  to  *hange  the  ilefinitimis  of  the  llnx 
functimis  from  those  given  in  (2)  to  thos*’  given  in  (12).  .Not*’  that  the  sidution  for  th*'  state 
is  tieede*!  in  order  to  evaluate  tlie  flux  functions  of  (12). 

•Next,  note  that  (11)  an*l  (2d)  are  identical  tlifferential  (^(ptations.  Thus,  the  chang*‘s 
tnatle  t*j  the  code  iti  or*l*‘r  to  treat  (11)  can  also  l>e  used  to  tr*’at  (2.1).  In  fact,  as  lotig  as  the 
ililferential  e(|uati«»n  ami  any  **ther  part  of  the  |>robl*‘m  sp*’cificaticui  ilo  not  exj^lir  itiy  tlep*'ntl 
on  th*’  ih’sign  |jaramet*’rs.  th<‘  analogous  r<’lation.s  will  b*’  th*’  sam**  for  all  th*‘  sen.-^it  i vilies. 

The  only  chang*’s  that  vary  fr*jm  one  sensitivity  calculation  to  anotlu’r  ar*’  thos*-  that 
arise  fr*un  conditi*uis  in  which  the  design  param*’t*’rs  app*’ar  *’Xplit  itly.  In  **ur  *‘xam|de.  for 
th*’  sf’iisitivity  with  r«’spect  to  .l/f,.  one  must  change  the  portion  of  th*’  co*!**  that  tn-afs  the 
inflow  ron*liti*jns  (  })-(.'))  so  that  it  can  iti.stead  treat  (1X)-(19).  In  th*’  probh’in  *  (jiisitU’n’il 
here,  the  nature  (i.e,  what  variables  are  specified)  of  the  b*>un*lary  *'onditions  at  th*’  inflow, 
ami  everywhere  «’lse.  is  not  affecte*!.  .N*>te  that  for  th*’  sensitivity  with  r<’spect  to  th*’ 
b*nin*lary  comliti<ui  (i-'j)  on  the  forebody  is  th*’  same  as  that  for  th*’  stat*'.  given  by  (b). 

For  the  ,sensitiviti*’s  with  respect  to  the  foreb(r«ly  <l*’sign  paramet*’rs.  th*’  inflow  boiimlary 
comlititnis  simplify  to  (27).  i.»’.,  tln’V  b*’c*un*’  homog*’n*'«ius.  Th**  boumlary  <(jmlition  at  the 
foreb<Mly  is  now  given  by  (dl)  or  (d-d).  Once  again,  the  nature  *T  th*’  boun*lary  lomlitions 
is  unrhange*!  from  that  for  the  state  ami  only  the  sp*’citi*’fl  data  is  difh’rent.  For  the  inthnv 
b(nin*lary  romlitif>ns.  we  may  still  specify  the  same  c*)mliti<jns  f(jr  the  s*’nsitivities.  but  now 
they  would  be  homogenetnis.  The  boundary  r*ui*liti*)ns  along  th*’  ffjrebody  chang*'  in  that 
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tlu'v  htH'otiu'  iiihoniog(Miet)us,  (rompare  (29)  and  (35)). 

In  smnniarv.  to  change  a  code  for  the  state  so  that  it  also  handles  the  sensitivities,  one 
must  rt'dtdine  the  tlnx  functions  in  the  differential  eciuations,  and  the  data  in  the  boundary 
conditions,  'fhe  changes  necessary  in  the  code  to  account  for  any  particular  relation  that 
(hx's  not  explicitly  involve  the  design  parameters  are  independent  of  which  sensitivity  one  is 
presently  c(jnsidering. 

Tlu'  previous  remarks  are  concerned  only  with  the  changes  one  must  effect  in  a  state  code 
in  otxh'r  to  handle  the  fact  that  one  is  discretizing  a  different  problem  when  one  considers  the 
sensitivities.  We  have  .seen  that  the.se  changes  are  not  major  in  nature.  However,  there  are 
additional  changes  that  may  be  needed  when  one  attempts  to  solve  the  discrete  eciuations. 
In  the  nunu'rical  n'sults  pre.sentc'd  below  we  use  the  finite  difference  code  “PARX!”  (see  [4] 
and  [N])  to  solv<>  the  state  and  .sensitivity  ecpiations.  However,  the  following  comments  apply 
eoually  well  to  other  X'FL)  codes  of  this  type. 

.Since  we  are  interested  in  steady  design  problems,  the  time  derivative  in  (1)  is  considered 
only  to  provide  a  means  for  marching  to  a  steady  state.  Now,  suppose  that  at  any  stage  of  a 
(iauss-.Nc'wton.  or  otiu'r  iteration,  we  have  used  PARC  to  find  an  approximate  steady  state 
solution  of  ( 1 )  plus  boundary  conditions.  In  order  to  do  this,  one  has  to  solve  a  .sequence  of 
liiK'ar  alg('braic  systems  of  the*  type 


(/  +  A/d(X^l:‘>))  =  (g};*’  +  AtB(QP))  ,  n  =  0, 1,2, . . . ,  (41) 


where  the  s<'(|U<Mice  is  terminated  when  one  is  satisfied  that  a  steady  state  has  been  reached 
and  where  deuotc's  the  discrc’te  approximation  to  the  state  Q  at  the  time  t  =  nAt.  We 
denote  this  stc-ady  state  solution  for  the  approximation  to  the  state  by  Q/,.  (4ne  problem  of 
the  type  (41)  is  solved  for  every  time  step.  In  (41),  the  matrix  A  and  vector  B  arise  from 
the  s])atial  discretization  of  tlu'  fluxes  and  the  boundary  conditions.  Both  of  tlu'se  depend 
on  the  state  at  the  prc'vious  time  level. 

Having  com|)Uled  a  steady  state  solution  by  (41),  the  task  at  hand  is  now  to  compute  the 
sensitivities.  We  will  focus  on  Q',  the  sensitivity  with  respect  to  the  inflow  Mach  number. 
.'\nalogous  results  hold  for  the  sensitivities  with  respect  to  the  forebody  design  parameters. 
Rec  all  that  given  a  state,  the  sc'usitivity  c'cpiations  are  linear  in  the  sensitivitic's.  Therefore*, 
if  OIK*  is  intc'rested  in  the*  steady  state  sensitivities,  instead  of  (11)  om*  may  directly  treat  its 
stationarv  version 


,Fl  ^  ^ 


(4'>) 


().r  Oy 

Since  (42)  is  linear  in  the*  compoiK'iits  ol  Q' ,  om*  dot*s  not  need  to  consi<h*r  marching  algo¬ 
rithms  in  order  to  compute  a  steady  s<*nsitivity.  ()m*  m(*rely  discr<*tiz«'s  (42)  and  solves  the 
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resultant  linear  system,  which  has  the  form 

A{Qu)Qu  =  (43) 

where  Q'/^  denotes  the  discrete  approximation  to  the  steady  sensitivity.  The  matrix  A  and 
vector  B  differ  from  the  A  and  B  of  (41)  because  we  have  discretized  different  differential 
equations  and  boundary  conditions.  Note  that  A  and  B  in  (43)  depend  only  on  the  steady 
state  Qi^  and  thus  (43)  is  a  linear  system  of  algebraic  equations  for  the  discrete  sensitivity 

Ql 

The  cost  of  finding  a  solution  of  (43)  is  similar  to  that  for  finding  the  solution  of  (41)  for  a 
single  value  of  n,  i.e.  for  a  single  time  step.  The  differences  in  the  assembly  of  the  coefficient 
matrices  and  right-hand-sides  of  (41)  and  (43)  are  minor.  Thus,  in  theory  at  least,  one  can 
obtain  a  steady  sensitivity  in  the  same  computer  time  it  takes  to  perform  one  time  step  in  a 
state  calculation.  If  one  wants  to  obtain  all  the  sensitivities,  e.g.,  K  -|-  1  in  our  example,  one 
can  do  so  at  a  cost  similar  to  ,  e.g.,  K-f-1  time  steps  of  the  state  calculation.  This  is  very 
cheap  compared  to  the  multiple  state  calculations  necessary  in  order  to  compute  sensitivities 
through  the  use  of  difference  quotients. 

Although  (43)  is  in  theory  no  more  complex  than  one  time  step  in  (41),  we  can  solve 
(42)  by  using  the  same  iterative  (or  another)  scheme.  The  simplest  approach  (but  certainly 
not  the  optimal  approach)  is  to  use  the  PARC  code  to  solve  (42)  by  time  marching.  In 
particular,  assume  that  Q/”*  is  a  solution  to  (41),  then  the  system 

[/  +  AM'CO!,"’)!  («')!:*■"'’  =  ((O')!"’  +  AiSWl,"')]  (44) 

can  be  used  to  find  given  ■  Thus,  one  makes  an  initial  guess  for  and 

and  then  iterates  (41)  and  (44)  simultaneously.  Also,  the  same  scheme  can  be  used 
to  compute  any  Qi-  =  i.e., 

[/  +  am'(c!;‘')]  (c3t)r"  =  [(<?»)!:'  +  a(f'(v!;‘>)]  ,  (45) 

In  practice,  these  “optimal”  estimates  of  speed  up  are  rarely  achieved.  Moreover,  as 
noted  above,  it  is  important  to  note  that  finite  difference  (FD)  and  sensitivity  equation  (SE) 
methods  do  not  necessarily  produce  the  same  results.  Since  the  ultimate  goal  is  to  find  useful 
and  cheap  gradients  for  optimization,  the  most  important  issue  is  whether  or  not  the  SE 
method  combined  with  an  optimization  algorithm  produces  a  convergent  optimal  design  as 
fast  as  possible.  We  have  tested  this  scheme  on  the  forebody  design  problem  and  the  next 
section  contains  a  summary  of  these  results. 
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6  An  Optimal  Design  Example 

In  order  to  illustrate  the  SE  method  and  to  test  its  use  in  an  optimization  problem,  we  used 
the  PARC  code  as  described  above  to  compute  sensitivities  and  the  used  these  sensitivities 
in  a  BFCIS/Trust  Region  scheme  to  find  an  optimal  shortened  forebody  simulator.  As  shown 
in  Figure  2.2,  data  was  generated  by  solving  the  Euler  equations  over  the  long  forebody  at 
a  Mach  number  of  2.0.  The  objective  is  to  find  a  forebody  simulator  with  length  one  half  of 
the  long  forebody  and  such  that  the  resulting  flow  matches  the  data  as  well  as  possible,  i.e. 
minimizes  J  along  the  outflow  boundary. 

The  shortened  forebody  was  parameterized  by  a  Bezier  curve  using  two  parameters. 
Thus,  there  are  three  design  parameters  p  —  {Mq,  P\,  Pi)-  The  algorithm  used  in  this 
numerical  experiment  was  based  on  using  the  PARC'  code  to  simultaneously  march  to  the 
steady  state  solutions  of  the  flow  and  sensitivity  equations.  We  made  no  attempt  to  optimize 
the  algorithm  since  the  main  goal  was  to  test  for  convergence. 

The  design  algorithm  proceeds  as  follows.  First,  an  initial  guess  for  the  optimal  design  is 
made,  i.e.,  we  select  a  A  good  selection  of  initial  parameters  can  be 

made  knowing  the  operating  conditions  of  the  aircraft  and  some  rough  guess  of  the  shape 
from  the  aircraft  forebody.  In  our  example,  we  chose  as  the  inlet  Mach  number  from 
the  computation  which  generated  our  data.  The  initial  guess  for  the  parameters  were  those 
used  to  generate  the  long  forebody  (although  corresponding  to  different  x-locations).  These 
parameters,  p°,  are  used  to  generate  a  grid,  the  inflow  and  forebody  boundary  conditions 
for  both  the  flow  (1)  and  sensitivity  equations  ((II)  and  (23))  and  an  initial  guess  for  both 
and  '•  our  example,  a  rough  guess  for  the  flow  field  uses  the  constant 

inflow  boundary  condition  throughout  the  flow  domain.  Likewise,  the  initial  guess  for 
is  taken  as  the  inflow  boundary  conditions  (given  in  equation  (18))  throughout  the  flow 
domain.  The  initial  guess  for  (Qa,)/”^  is  initially  taken  as  zero  (except  on  the  forebody).  The 
systems  (41),  (44)  and  (45)  are  then  solved  simultaneously  (in  our  case  the  left  hand  side 
matrix  is  the  same  for  (41)  as  for  the  sensitivity  equations  (44)  and  (45),  i.e.  A  =  A')  for 
the  updated  \  ^  (^)/  updated  is  then  used  to  formulate 

(41),  (44)  and  (45)  and  .solve  for  (Q/i)^”"^*^  and  ^  Then  one  iterates  until  the 

desired  convergence  is  achieved.  In  our  example,  the  residuals,  AQ/,  =  were 

converged  to  approximately  10“’®  (in  800  time  steps).  The  outflow  data  Q/^  and  3.re 

then  used  to  compute  J{p^)  and  Vj7(p‘’). 

The  optimization  algorithm  consisted  of  a  BF(«S  secant  method  coupled  with  a  “hook” 
step  model  trust  region  method  [5].  The  initial  Hessian  was  obtained  by  finite  differences 
on  VJ{p).  The  function  and  gradient  information  needed  by  the  optimization  algorithm  is 
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obtained  by  calling  the  modified  PAR(!  code  with  p  =  p. 

This  algorithm  was  tested  for  the  case  where  the  forebody  simulator  was  allowed  to  have 
the  full  length  of  the  body  generating  the  data.  In  this  case  the  optimization  algorithm 
produced  exact  data  fits,  i.e.  J[p’')  —  0  and  it  recovered  the  parameters  used  to  generate 
the  data.  However,  the  more  realistic  test  (constraining  the  length  of  the  forebody  simulator) 
also  produced  a  convergent  design  and  reduced  the  cost  functional  significantly. 

Figure  6.1  shows  the  flow  field  over  the  long  forebody.  Observe,  that  there  is  a  shock  in 
the  flow.  As  noted  in  [2],  shocks  can  cause  difficulties  if  one  is  not  careful  in  the  .selection  of 
an  appropriate  numerical  scheme.  High  order  schemes  can  produce  (numerically  generated) 
local  minimum  that  can  cause  the  optimization  loop  to  fail.  This  problem  is  avoided  here 
because  the  numerical  viscosity  in  PARC  (required  for  stability)  is  sufficient  to  “smooth” 
the  cost  functional  (see  [2]  for  details). 

Figure  6.2  shows  the  shape  and  flow  field  of  the  optimal  shortened  forebody.  This  design 
was  obtained  after  12  iterations  of  the  optimization  loop.  Figures  6. 3-6. 6  show  the  1*‘,  2’“^, 
3”'^,  and  12‘^‘  iterations  for  each  of  the  flow  variables.  The  initial  guess  for  the  parameters 
were 

/=  ((a/o')°,P,,P,)  =(2.0,0.10,0.15) 

and 

J(p“)  =  3.2339. 

The  “converged”  optimal  parameters  are 

p-  =  =  (2.020,0.294,0.156) 

with 

J(p-)  =  0.2229. 

Observe  that  the  cost  function  was  decreased  by  more  than  93%.  Figures  6.7-6.10  show  a 
comparison  of  the  flow  fields  for  the  optimal  shortened  forebody  simulator  and  the  data.  The 
optimization  loops  converged  rapidly.  For  example,  J[p^)  —  0.2334  and  J{p^)  =  0.2289. 
This  is  due  to  the  fact  that  the  shock  location  was  found  quickly. 

Note  that  although  the  flows  are  close,  there  is  a  significant  error  near  the  forebody.  This 
can  also  be  seen  in  the  plots  in  Figures  6.11-6.14.  It  is  worthwhile  to  note  that  the  match 
is  good  considering  the  fact  the  shortened  forebody  is  constrained  to  be  one  half  the  length 
of  the  “real”  forebody  and  only  two  Bezier  parameters  are  used  to  model  r(  ).  It  is  also 
important  to  note  that  the  shock  is  captured  by  the  optimal  design.  In  particulai  serve 
in  Figures  6. 3-6. 6  how  the  optimization  algorithm  “shapes”  the  shortened  forebod  >  that 
the  optimal  shape  has  a  blunt  nose.  This  is  necessary  in  order  to  generate  the  correct  shock 
location  at  the  outflow. 
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7  Conclusions 


The  numerical  experiment  above  illustrates  that  the  SE  method  can  produce  sensitivities 
siiitable  for  optimization  based  design.  There  are  a  number  of  interesting  theoretical  issues 
that  need  to  be  addressed  in  order  to  analyze  the  convergence  of  this  approach.  Moreove* 
one  should  investigate  “fast  solvers”  for  the  sensitivity  equations  (multi-grid,  etc.)  cis  we! 
develop  numerical  schemes  that  are  not  only  fast,  but  produces  consistent  derivatives  wh- 
possible. 

Finally,  we  note  that  we  have  conducted  a  number  of  timing  tests  which  compute  sen¬ 
sitivities  to  compare  the  SE  method  with  the  finite  difference  method.  In  particular,  we 
observed  that  for  the  problem  above  (with  three  design  parameters),  the  SE  method  needed 
only  58%  of  the  CPU  time  required  by  finite  differencing.  When  twenty  design  parameters 
were  used,  the  SE  method  produced  these  sensitivities  in  about  38%  of  the  time  required  by 
finite  differencing.  These  early  numerical  results  indicate  that  considerable  computational 
savings  may  be  possible  if  one  extends  and  refines  the  basic  SE  method  presented  here. 
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Figure  2.1: 2D  Forebody  Problem 
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Figure  2.2:  A  2D  Optimal  Forebody  Design  Prob!  ;  ; 


Figure  6.1:  Long  Forebody  (Outflow  Data  to  be  Matched) 
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Figure  6.3:  Iteration  to  Optimal  Forebody  Design: 
Density 
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Figure  6.5:  Iteration  to  Optimal  Forebody  Design: 
X-Component  of  Momentum 
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Figure  6.6:  Iteration  to  Optimal  Forebody  Design: 
Y-Component  o£  Momentum 
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Figure  6.7:  Comparison  o£  Optimal  Shortened  and 
Long  Forebody:  Density 
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Figure  6.8:  Comparison  of  Optimal  Shortened  and 
Long  Forebody:  Energy 
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Figure  6.9:  Comparison  of  Optimal  Shortened  and 


Long  Forebody:  X-Component  of  Momentum 
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